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Abstract. We propose a method for the fast generation of a quantum register 
of addressable qubits consisting of ultracold atoms stored in an optical lattice. 
Starting with a half filled lattice we remove every second lattice barrier 
by adiabatically switching on a superlattice potential which leads to a long 
wavelength lattice in the Mott insulator state with unit filling. The larger 
periodicity of the resulting lattice could make individual addressing of the 
atoms via an external laser feasible. We develop a Bose-Hubbard-like model 
for describing the dynamics of cold atoms in a lattice when doubling the 
lattice periodicity via the addition of a superlattice potential. The dynamics 
of the transition from a half filled to a commensuratcly filled lattice is analyzed 
numerically with the help of the Time Evolving Block Decimation algorithm and 
analytically using the Kibble-Zurek theory. We show that the time scale for the 
whole process, i.e. creating the half filled lattice and subsequent doubling of the 
lattice periodicity, is significantly faster than adiabatic direct quantum freezing of 
a superfluid into a Mott insulator for large lattice periods. Our method therefore 
provides a high fidelity quantum register of addressable qubits on a fast time scale. 
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Figure 1. (a) The initial profile of the lattice is associated with the value of 
s = 1. The periodicity of the lattice is then progressively doubled (b) until it 
reaches its final profile (c) associated with the parameter s = 0. The number 
of sites M corresponds to the number of unit cells in the large lattice limit. By 
starting with a filling factor of n = 1/2, this procedure leads to a lattice with 
filling factor of n = 1. 



1. Introduction 

Systems of cold atoms trapped in optical lattices provide the unique opportunity to 
coherently manipulate a large number of atoms [TJ [3] . The remarkable degree of 
experimental control offered by these systems, as well as the possibility to use the 
internal hyperfinc states of the atoms to encode qubits, make them particularly suited 
for quantum information processing (QIP). In this context, optical lattices in a Mott- 
insulating (MI) state with unit filling can be viewed as the realization of a quantum 
register, and it is possible to collectively manipulate the qubits stored in such a register 
experimentally [4l|5]. However, in many quantum computing schemes based on neutral 
atoms stored in optical lattices the application of single qubit gates [3] or single qubit 
measurements [6] requires the ability to address single atoms with a focused laser 
beam. This remains experimentally challenging since these operations have to be 
performed without perturbing the state of other atoms in their vicinity. 

A number of strategies have been proposed to circumvent this problem by using 
global operations [ZlIHlIS], e.g. via "marker atoms" which are moved to a particular 
lattice site and interact with the corresponding register qubit such that an external 
laser affects only that qubit [lOj . Another way is simply to use a quantum register 
in which the atoms are distant enough such that they can be addressed individually 
by a laser. This method requires an optical lattice with filling factor n = 1 in the MI 
state with a sufficiently large distance between the atoms [IT] . The initialization time 
of a MI state is proportional to the tunneling time of the atoms between neighboring 
sites. Therefore, by using the conventional method of quantum-freezing a superfluid 
(SF) state [Hiniinilll], it scales exponentially with the lattice spacing [T5 l [TB I fT7] . 

In this paper we propose an alternative method to generate a long wavelength 
lattice with one atom per site in the MI state. Starting with a one dimensional lattice 
with a short period — and hence a short initialization time — and filling factor n = 1/2 
we remove every second potential barrier by adiabatically turning on a superlattice. 
This superlattice potential has already been experimentally realized [HI [19]. This 
procedure eventually leads to a long wavelength lattice in which the periodicity has 
been doubled and where the atoms are in a MI state with n = 1 (see figure [T]). 
This scheme does not require changing the angles of the intersecting laser beams. 
Furthermore, we show that our method allows for the initialization of a MI state 
with unit filling factor on a time-scale which, although scaling exponentially with the 
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final lattice spacing, is approximately one order of magnitude smaller than the direct 
quantum-freezing method. Although we only consider the case of an homogeneous 
lattice, the results presented in this paper extend to the case of weak harmonic 
confinements, quartic [50] and box traps pT| . 

This paper is structured as follows: In Sec. [2] we introduce the model used 
to describe the system dynamics. In Sec. 13.11 we discuss ground state properties, 
particularly two-site correlation functions and quasi momentum distributions. Also, 
we present and discuss numerical results for the probability of staying in the ground 
state during the transition depending on the speed of the ramping. In Sec. 13.21 we 
apply the analytical Kibble-Zurek theory and compare it with our numerical results. 
Finally, we summarize and conclude in Sec. ID 

2. Model 

We consider a gas of interacting ultracold bosonic atoms loaded into a three 
dimensional optical lattice. The lattice is formed by pairwise orthogonal standing 
wave laser fields and its optical potential is given by |22j 

Vodr, s) = Vsix, s) + Vt [sin^{kz) + sin\ky)]. (1) 

Here Vt is the depth of the potential in the y~ and z-directions created by pairs of 
lasers with wave number k = 27r/A, wave length A and period = A/2. The extension 
to the case of different optical potentials in the y- and z-directions is straightforward. 
In the x-direction two pairs of laser beams with a long wavelength A^ and short 
wavelength As = Al/2 are applied. The potential in the x-direction is thus given by 

Vs{x, s) = (1 - s) sin^ (fcix) + s sin^ (fc^x) , (2) 

with Vq the depth of the lattice, k^ = 2tt/Xl and ks = 2tt/Xs. The depths of the 
potentials will be expressed in units of the recoil energy Efj = k\/2m with m the 
mass of the atoms (taking h—1 throughout). The parameter s G [0, 1] is determined 
by the relative intensities of the two pairs of lasers. By changing s from 1 to the 
lattice profile is continuously transformed from a sinusoidal potential with a small 
period as = Xs/2 to one with a long period a = Ai/2, thus halving the number of 
lattice sites per unit length along the x-direction (see figure [T]) . The lattice constant 
a corresponds to the size of a unit cell for s < 1 [J. We refer to the lattice profile with 
parameters s = 1 and s = as to the small lattice limit and the large lattice limits 
respectively. 

The Hamiltonian of the system in second quantization reads 

y"dr^t(i.)/jQ(f)v3(r) + | y"dr7At(r)^t(i.)^(r)^(r), (3) 

where /io(f) = -(l/2m)V2 + VbL(f) is the one-particle Hamiltonian. The symbol 
f = (r, s) represents the position variable r and lattice parameter s. The interaction 
between the atoms is modelled by s-wave scattering with g = Anas/m where is the 
s-wave scattering length. The bosonic field operators obey the usual commutation 
relations [tp{r),il'^ {r')] = S{r — r') with 6 denoting the Dirac delta function. 

We restrict our considerations to the case where Vt is sufficiently large so that 
motion along the y- and z-directions is frozen. The dynamics of the system is then 
effectively one dimensional along the x-direction. As shown in figure[2^ the two lowest 

f To avoid any discontinuity, we work with a lattice periodicity of a even in the case s = 1. 
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Figure 2. Band structure along the direction in (a) the small and (b) 
large lattice limit for Vb = lO-E/j. The points represents the values of q for 
q = 0, 7r/a, 7r/2a. The value of Sa,q correspond to the energy of single-particles 
with momentum q in the a— th Bloch band. In the small lattice limit, the two first 
Bloch bands arc connected. 



bands of the Hamiltonian /iq.x = — (l/2m) {A/ Ax) + Vs(a;, s) are separated in energy 
by much less than the typical motional excitation energy iJex = V^VqEji for 
values s ~ 1 Therefore, despite assuming that atoms loaded into the lattice are 
ultracold, we have to consider the two lowest Bloch bands in x-direction to obtain an 
accurate description of the atomic dynamics. However, excitations to higher bands in 
the y- and z-directions are neglected in our investigations since we assume that the 
temperature of the atomic cloud is ksT <C -Ecx- We then expand the bosonic field 
operator as 

M M 

ip{r) = ^ (t)a4r) a, +^ 4>b4^)k. (4) 

1=1 i=l 

where M is the number of lattice sites and a| {a — a, b) creates a particle in the 
mode associated with the localized function (f>a^i{r) centered at site i. The mode 
functions (/)Q_i(f) = Wa,iix, s)Wi,o{y)Wi,o{z) are factorized into a product of well 
localized Wannier functions (WF) Wi,o of the lowest Bloch band in the y- and z- 
directions [23j and mode functions Wa,i{x,s) in a;-direction. 

The aim of the next section is to describe the single particle dynamics in the tight 
binding (TB) approximation. If we were to use Wannier functions for Wa^i{x,s) this 
approximation would restrict our model to sinusoidal Bloch bands [24] . Because of the 
deviation of the lowest two bands from a sinusoidal dispersion relation (see figure [5^) 
when s w 1 we instead use generalized Wannier function (GWFs) for Wa^i{x,s) (see 
[Appendix A| for a detailed definition and a description of their properties) [25| . By 
exploiting the optimization procedure described in [Appendix A[ we calculate GWFs 
Wa^i{x,s) which are well localized at lattice sites i. Typical shapes of GWFs and 
the effects of optimizing their localization are shown in figure [3| We note that these 
GWFs are in general composed of superpositions of Bloch orbitals of both bands and 
are not related to Wannier functions by a local transformation. Only when s ~ 
the Wa,i{x,s) are equivalent to the Wannier functions of the first and the second 
Bloch band, respectively. Finally, they are always (anti-)symmetric with respect to 
the center of the lattice site i for a = a {a = b). 

§ The two lowest bands form two segments of the lowest Bloch band if a cell size of a/2 is used in 
the case s = 1. 
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Figure 3. (a) Optimized (thick line) and non-optimized (dashed hne) GWFs 
associated with the first mode for Vq = SO-E/j and the lattice profiles s = 1. The 
area between the optimized and non-optimized mode functions is shaded in order 
to illustrate how the localization procedure reduces the spread of the optimized 
function, (b) The square of the mode functions for Vq = SOEn and s = 1. By 
combining the mode functions shown in (b), we can construct two new mode 
functions corresponding to particles localized in either the left (c) or the right (d) 
well of a given site. 



2.1. Single-particle Hamiltonian 

Inserting the approximate field operator equation (|4|) into the first term of equation ([3]) 
yields the single-particle part of the Hamiltonian in terms of a] and hi . Applying the 
tight binding approximation, which amounts to keeping only nearest-neighbor hopping 
terms, the single-particle Hamiltonian can be approximated by 

Hois) = ^ (^Jfcb(s)5j6i+i - Jaais)alat+i + h.cj 



=1 



i^Jbais) bloi+i - Jab{s) albi+i + h.c. 

1=1 

M 

Y,(Va{s)ala, + n{s)blk), (5) 



where 



Japis) = J dxw'^ i{x,s)ho,x{s)wf3^i+i{x,s), (6) 

is the hopping matrix element between neighboring sites along the x-axis and 

Va{s) = J dxw*^ i{x,s)ho^x{s)wa,i{x,s), (7) 

is the local on-site energy of a particle in mode a. Note that hopping between modes 
a and b within one site is not allowed by the symmetry properties of the GWFs. 
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Figure 4. Parameters of the effective Hamiltonian -ffcff as functions of s for a 
lattice depth of Vq = SOEji and Vt = QOEji with d = a o3p . The hopping matrix 
elements (a) have been calculated using the method described in [Appendix B[ 
while the on-site interaction energies (b) are calculated using optimized GWFs. 



However, the inclusion of non-zero hopping matrix elements Jab and J^a is essential 
to accurately reproduce the single particle behaviour of the full Hamiltonian ([3]). The 
symmetry properties of WFs would not allow the inclusion of these terms [53] . 

For periodic boundary conditions, the parameters V„ and Ja^p can be found 
from the band structure without explicit calculation of the mode functions (see 
[Appendix B[ ). The numerical values of Va and Ja^p for Vq = SOi^i^ are shown in 
figure 3^. Using these parameters, we find that Ho{s) very accurately reproduces 
the band structure of the exact Hamiltonian for all values of s, thus justifying the 
utilization of the TB approximation and the corresponding GWFs. 



2.2. Interaction Hamiltonian 

To calculate the interaction matrix elements of H the explicit form of the localized 
GWFs is needed. We find that for Vb > IOEr, off-site interaction terms are at least 
two orders of magnitude smaller than on-site interactions. We therefore only keep the 
dominant on-site terms and find the interaction Hamiltonian Hj 

Hjis) ^j:^^n^ {nl _ 1) + Mf) n\ (n\ - l) 



1=1 



+ E ^y^ + iMa^d, + dld^J^k) , (8) 

i=l 

where n° = aja, and = b\bi are the site-occupation number operators. The on-site 
interaction matrix elements are given by 



UaAs)^9 j dvwl i{r)w*p^^{T)wc,i{v)wp^,{r). (9) 

The numerical values of Ua^p as a function of the lattice profile s are shown in figure[3)3. 
Figure [41d shows that their values become equal as s ^ 1 for sufficiently large Vq. 

Combining the single- and two-particle contributions the effective Hamiltonian 
describing the system dynamics is given by 

Heeis)=Hois)+Hiis). (10) 

By using Hcff{s) for s varying in time we implicitly assume that the system 
adiabatically follows changes in the mode functions Wa^i- For all dynamical 
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calculations carried out in this work we have carefully chosen the time dependence 
of s so that such non-adiabatic contributions can safely be neglected. 



2.3. Limiting cases 

In the small lattice limit (s = 1), the superpositions ipL,i — {wa,i — wi,_i)/\/2 and 
V-R.i = {wa.i + Wb.i)/\/2 correspond to mode functions localized in the left and in the 
right well of site i respectively (see figure[31; and[3Jl). The associated bosonic operators 
are defined by 

ch = ^{a\-h% ct,,, = i=(at + 6l). (11) 

Given that the new mode functions are sufficiently localized within each well, the 
parameters of H^s for s = 1 can be written as 

Va=E-J, Vb^E + J, U^,p = U/2, (12) 

where J = / dx lp\ .Jio,x^la+i-, E ^ J da; (/?2,i^o,a:'^L,i and U = g J dx \ipB^A'^- Notice 
that the parameters shown in figure [3] are consistent with these equations. Expressing 
the Hamiltonian (jlOp in terms of the operators cj^ ^ (a = L, R) and using the 
parameters 1121 we find that 

2A/-1 2M jj 2M 

Ks = II + h.c. +eY, 44c»'C,', (13) 

i' — 1 i' — 1 i'—l 



where 



if i odd, 
if i even. 



4 = <! .r^ .„ . (14) 



Therefore, as expected, Hes{s = 1) is equivalent to the one-band Bose-Hubbard 
model (BHM) ^ with 2M sites. 

The mode functions keep their symmetry with the two peaks moving towards 
the center of the cell when s is decreased. When s « is reached Wa,i{x,s) and 
Wb^i^Xjs) are equivalent to the Wannier functions of the first and the second Bloch 
band, respectively. In this limit we obtain a standard two band Bose-Hubbard model 
for M sites. 



3. Time-scale for the preparation of the quantum register 

In this section we present numerical as well as analytical results characterizing the 
ground state properties of the system and the time-scale necessary to initialize the 
quantum register. 



3.1. Numerical results 

The numerical calculations have been carried out using both the exact matrix 
representation of the Hamiltonian iJofr and the Time-Evolving Block Decimation 
(TEBD) algorithm (see [Appendix C| for details). 

We have evaluated the ground state and dynamical properties of our system for 
two different values of g corresponding to different interaction regimes. These values 
have been chosen such that in the large lattice limit (with filling factor 7i — 1) the 
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ground state is a Mott-insulator state for gi and g2 with Jaa/Uaa < 0.3 [^Hj. In 
the small lattice limit, gi and g2 produce ground states corresponding to a strongly 
interacting Tonks- Girardeau (TG) gas {U/J = 215) and a superfluid {U/J = 7), 
respectively. 



3.1.1. Ground states properties In the large lattice limit, the ground state IV'o.l) of 
the system is populated exclusively by particles in the lowest Bloch band, i.e. a-mode 
particles. Therefore, in this limit, we only consider the one-particle density matrix 
given by 

= (V'|ala,#), (15) 
where \'>p) is the state of the system. The quasi- momentum distribution of particles in 
the a-mode is given by [17] 

M 



W = E ^-"""^'^'^ pfM)- (16) 



As expected, in this limit and for commensurate filling n = 1 the ground state 
is a Mott-insulator for both values of g (see figures [S^-d). The quasi-momentum 
distributions show that in this limit particles arc uniformly distributed in the first 
band (see figures [5^ and [5]:). 

In the small lattice limit, the one-particle density matrix is given by 

pf,^,i^j) = {^P\4,c,,\i'), (17) 

where the operators c]^, are constructed from the a] and hi operators using the 
transformations pi[l . The quasi-momentum distribution is given by 

2M 

= ^ E e-'"^^'-''^ pf'fW- (18) 

In this limit (with filling factor n — 1/2), the characteristics of the ground state 
IV'o.s) depend on the value of g. For g = gi, the ratio U/J~ 215 and the ground 
state's correlations as well as the quasi-momentum distributions (see figures [l^-f) are 
characteristic of a TG gas (see e.g. [28] and references therein). The value g = 52 
yields the ratio U/J = 7 and the system exhibits the behaviour of a superfluid (see 
figures [5^-h) , with particles occupying mainly the q = momentum state. Notice 
that one-dimensional systems described by the BHM with a fixed filling factor n = 1 
cross the MI-SF phase boundary at the critical point {J/U)c ~ 0.3 [2^. Thus, in 
our case the superfiuid behavior of the system in the small lattice limit is due to the 
fractional filling of the lattice. 



3.1.2. Simulation of the dynamics Starting from a system with half-filling and a 
lattice profile s = 1, we investigate the time-scale required to obtain a nearly perfect 
MI state — or quantum register — with filling factor n = 1 by ramping the lattice profile 
down to s = 0. The quality of the register is determined by the fidelity 

F=M^o.l)\\ (19) 
defined as the overlap between the state of the system \tp) at the end of the ramping 
process and the ground state in the large lattice limit IV'o.l) • Furthermore, we calculate 
the fluctuations of the number of particles in the a-mode at site i 

An,,. = ^((^^iV^, (20) 
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Figure 5. Ground state one-particle density matrix and quasi-momentum 
distribution for M = 12 lattice sites and the parameters of figure |4] The figures 
(a— d) correspond to the large lattice limit and the figures (e— h) to the small lattice 
limit. The quasi-momentum distribution (a), (e) and one- particle density matrix 
(b), (f) are for g = g\. The quasi-momentum distribution (c), (g) and one- particle 
density matrix (d), (h) are for g = g2- 
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Figure 6. Dynamical simulation of -ffgff using the parameters shown in figure |4] 
for two diff'eront values of g. The ramping strategies Sgap and Sman are shown 
in figure [7^. Simulation results for the fidelity and particle-number fluctuations 
using: (a-b) a linear ramp and g = gi; (c-d) a linear ramp and g = 32; (e— f) 
the ramp Sgap and g = gi; (g— h) the ramp Sman and g = g2- The vertical lines 
indicate the value of the particle tunneling time (ttun ~ 1/Jaa) in the large lattice 
limit. 



Fast initialization of a high-fidelity quantum register using optical superlattices. 



11 



1 




^ 0.5 



inian 



D,2e 















t 



t 



Figure 7. (a) DifTcront ramps used in our numerical calculations, (b) Transition 
probability between the ground and the first (labeled by e) and second (labeled 
by 2e) excited states as a function of time for the ramps and Sman- The 
transition probabilities have been calculated via the exact diagonalization of Hgff 
for M = 4 and a quench time of tq = 20/ En. 



where (o) = (-01 o lip). Since particle-number fluctuations are suppressed in a MI state, 
non-zero fluctuations indicate the presence of excitations, such as double occupancies 
or particles in the second band, in the final state. 

We test different ramps by simulating the system dynamics between tj = to 
ti = '''Q where tq is the ramping time (the time required to complete the ramping 
process from s = 1 to s = 0). Here, each ramp a corresponds to a function s ~ Sa-{t) 
varying from Scr(O) = 1 to Sa{TQ) = 0. 

For linear ramping we use siin(i) = {tq — tj E-ii)Itq- The fidelity and particle- 
number fiuctuations obtained using this strategy for different quench times and g ~ g\ 
and g = gi are shown in figures [G^-d. The linear ramp is shown in figure [7^. 

Another ramping strategy we use consists of adapting the velocity of the ramp 
proportionally to the energy gap between the ground and the first excited state. This 
ramp is denoted by Sgap(i)- We evaluate the ramp function Sgap(i) numerically (see 
[Appendix D[ ) for a system with M = 4 sites and Vb = 30i?i^. The ramp Sgap(t) for 
g = g\ Ss shown in figure [7^. We expect this ramping strategy to be more efficient 
than the linear one, since accelerating the ramp when the gap is large while slowing it 
down when the gap is small should suppress transitions of particles to excited levels. 
Our numerical calculations have shown that when g = gi, the utilization of this ramp 
does indeed reduce the quench time needed to obtain a nearly perfect fidelity to a half 
of the tunneling time in the large lattice limit (see figure [6^-f). For 5 = 52, we find 
that compared to siin(i), this strategy only marginally improves the fidelity. 

In order to reduce the time required to obtain a given fidelity for g = g2, & more 
sophisticated ramping strategy is needed. In the following, we provide a simple method 
to estimate the efficiency of different ramps without running a complete dynamical 
simulation of the system. The transition probability between the ground and some 
excited state at time t for a ramp a is approximately given by |29j 





PSk{t)^A- (Ol^^lfc) [1 - cos(^o;.t)] , 



(21) 
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where \k) = \k{cr,t)) is the fc-th instantaneous cigcnstat^oi Hcft{sa-(t)) and ojofe is the 
transition frequency between the ground state and \k) . Therefore, the assessment of 
the efficiency of a ramp can be done by evaluating the functional 

1 , f^Q 

A{<j,tq)^—Y, dtPS,{t), (22) 

where the index k runs over all the values associated with levels connected to 
the ground state. The functional A{a,TQ) corresponds to the average transition 
probability per unit time for a given strategy a and quench time tq . We calculate the 
value of the functional A{a, tq) numerically via exact diagonalization of for a small 
system. This method allows to optimize ramps by minimizing the value of A{<7,tq). 
The optimized ramp for a small system is then used in the simulation of larger systems. 
For instance, the strategy Sman(i) shown in figure [7^ was designed and optimized 
manually using this method. For g = 92, we find that A(lin, tq) / A{ma'n, tq) w 2.3 for 
TQ = 20/ En (sec figure [Tj^), thus showing the better efficiency of the strategy Sman(i) 
compared to sii,i(t). As shown in figure [6l dynamical simulations of the system with 
g — g2 confirm that this strategy reduces the time required to obtain a given fidelity. 

For systems initially in the superfiuid regime {g ~ .92), the fidelity curves exhibit 
small oscillations (see figure [BJ2). These can be understood from time-dependent 
perturbation theory as oscillations occurring when some of the frequencies involved in 
the Fourier decomposition of the perturbation Hamiltonian enter into resonance with 
system frequencies. This is expected since supcrfluids have a dense spectrum at low 
energies and are therefore likely to enter into resonance with one of the frequencies 
of the perturbation Hamiltonian [30| . Hence, the amplitude of the oscillations in the 
fidelity curve associated with a ramping strategy si{t) should be smaller than those 
associated with a ramping strategy S2{t) if A{1,tq) < A(2,tq) for all tq. This is 
what is observed from our numerical simulations (sec figures |6f-h). 

The quasi-momentum distribution of the particles in the a-mode at the end of 
the different ramping processes for a system with g = gi are shown in figure [HI For the 
linear ramp, the quasi-momentum distribution of particles shown in figure[5^ does not 
correspond to that of a MI state for the quench times considered. In contrast, for the 
ramp Sgap the quasi-momentum distribution becomes approximately flat for quench 
times of tq > 200/ En. Even for the fastest ramps, we find that the occupation of the 
6-mode is less than 2% of the total number of particles. Thus, the experimental 
measure of the register fidelity can be made by comparing the quasi-momentum 
distribution of particles in the final state with, e.g. that shown in figure [5}:. 

3.1.3. Discussion of the numerical results In the BHM with U / J <^ 1, the tunneling 
time 1/J determines the adiabatic time-scale of the system. However, as soon as 
many-body interactions are sufficiently large, this time-scale often becomes very 
non-adiabatic [3T]. The main observation that can be drawn from the numerical 
calculations presented in the last section is that by preparing the system in a TG 
state {g = gi), and using an efficient ramping strategy, it is possible to initialize a 
very deep MI state on a time scale equivalent to half the tunneling time in the large 
lattice limit (see figure |6Jd) . The time required to initialize a MI state with unit filling 
as well as a TG state with half filling is approximately ten times the tunneling time of 

II Note that neglecting non-adiabatic changes of the GWFs does not correspond to P^i^it) = at all 
times. 
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Figure 8. The quasi-momcntum distributions at tlie end of a ramp for a system 
of M = 12 sites using the parameters shown in figure |4] for g = gi- (a) sii^; (b) 

Sgap • 

the final system [ini [SH [Ml US] ■ Since the tunneling time in the large lattice limit is 
two orders of magnitude larger than in the small lattice limit, the total time required 
to initialize a MI state using our procedure is an order of magnitude faster than the 
direct quantum freezing method. In this estimation we assume that the initial BEC 
has zero temperature, i.e., we do not take the effect of defects present in the initial 
state into account. 

The experimental realization of initial states with g = gi and g = g2 can be 
achieved using Fcshbach resonances. For a magnetic Feshbach resonance fluctuations 
in the magnetic field result in fluctuations of the size of the gap between the ground 
and the first excited state which will affect the performance of our scheme. For 
instance, in the case g = gi magnetic field fluctuations of 10 mG will change the 
gap by approximately 0.5% for ®^Rb or ^^'^Cs atoms [331 [M]. We assume adiabatic 
evolution of the system and thus these fluctuations will have negligible repercussions 
on the fidelity of the final state. To realize the superfluid regime with g = g2 very 
stable magnetic fields are required. Magnetic field fluctuations of e.g. 1 mG will lead 
to gap fluctuations of approximately 1% in ^"^Na and ®^Rb. We finally remark that 
our scheme could also be used without employing a Feshbach resonance. This case 
corresponds to an intermediate value of g between gi and 52 • While a detailed analysis 
of the intermediate regime is beyond the scope of the present work we do not expect 
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qualitative differences compared to the interaction strengths considered here. 
3.2. Analytical results 

In this section wc derive an approximate expression for the quench time required to 
obtain a given fidcUty in the case of a hncar ramp. 

The energy spectrum of systems in the TG and superfluid regime is gapless0(see 
e.g. [35j). while in the Mott-insulating regime, the gap between the ground and first 
excited state is proportional to the on-site interaction energy [30j . Since the relaxation 
time r(t) of the system — the time required by the system to adjust to a change of 
parameters at time t — is inversely proportional to the gap between the ground and 
the first excited state, the relaxation time in the small lattice limit is large, while it is 
small and finite in the large lattice limit. This observation suggests that the adiabatic- 
impulse (AI) assumption from Kibble-Zurek (KZ) theory can be used to evaluate the 
adiabaticity of a ramp with respect to the quench time |361 1371 138] . 

The AI approximation is based on the following considerations: (i) When the 
gap between the ground and the first excited state is large, the relaxation time of the 
system is short and thus a system starting its evolution in the ground state remains 
in the ground state, i.e. its evolution is adiabatic. (ii) When the gap between the 
ground and the first excited state is small, the system's relaxation time is large and 
the system no longer adapts to changes of the Hamiltonian's parameters and its state 
becomes effectively frozen. The system is then in the impulse regime. The instant t 
at which the system passes from the impulse to the adiabatic regime, and inversely, is 
defined by the equation [371 ISSl IM] 

r(t) = at, (23) 

where a = 0(1) is a constant. Note that t is a time and not an operator. In the AI 
approximation, the time-evolution of the system is either adiabatic or impulse. Thus, 
the density of defects '3 , which corresponds to the density of excitations caused by a 
change of parameters which drive the system from the impulse to the adiabatic regime 
can be approximated by [37j 

^C.|(*e(i)|*s(0))|', (24) 

where |^'g(0)) and |^'e(t)) arc the ground and first excited states at the initial time 
t = and at time i ~ t, respectively. Hence, without solving the time-dependent 
Schrodingcr equation it is possible to make predictions for the density of defects ([24]) 
resulting from a given dynamical process. 

In order to apply the KZ theory to our problem, wc develop an effective model 
describing the system dynamics. A similar model was recently examined by Cucchietti 
et al. [38] . We find from numerical calculations (see figure [7)d) that most of the 
excitations created in the system are caused by transitions from the ground to the 
first excited state. Furthermore, we examined the form of the eigenvectors of the 
Hamiltonian (jlOp in both limits for different system sizes via exact diagonalization. 
This revealed that both the ground state and the first accessible excited state can be 
approximated by an expansion in only two basis states. The elements of this reduced 
basis set are given by 

M 

|l)=(g)|l>., (25) 

i=l 

^ In finite size systems, tiie spectrum is not gapless, only very dense. 
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|2;0;1;---;1) + |2;1;0;1;---;1) + |2;1;1;0;1;---;1) 



+ \1;2;1;0;1;---;1) + ■■■)/ ^M{M-1), (26) 

where, e.g. |2; 0; 1; • • • ; 1) = |2)i«)|0)2®|l)3«)- ■ "^ll)*/ with |n), = (1/V^)(aj)"|vac). 
The basis state 1 2) corresponds to a superposition of all possible states of a system of 
M particles in the a-modes with (M — 2) singly occupied sites, one doubly occupied 
site, and one empty site. In the limit M —^ oo, the matrix representation Hfi of 
Hamiltonian (fTO]) in the reduced basis {|1), |2)} reads, up to a constant energy Va 



V2Jaait) Uaa{t) 

The instantaneous eigenstates of equation (^7)) associated with the energies of 
the ground and first excited levels are given by \g{t)) = — sin(^^(^)/2)|l) + 
cos(6'(t)/2)|2) and|e(t)) = cos(6'(t)/2)|l)+sin(6'(t)/2)|2), respectively, with cos(6'(t)) = 
—Uaa{t)/\Uaa{t)'^ + 8Jaa(i)^]^, G [0, tt]. Furthermore, we approximate the 
parameters Jaa and Uaa as linear functions of time 

Jaa{t) = ^Jaa{TQ - t) / Tq , 
Uaa{t)=Uinit+AUaat/TQ, 

where AJaa = |Jaa(s = 1) - ./qq(s = 0)1 , AUaa = |t/init " Uaa{s = 0)| and 

;7init = niin[C/aa]- 

We further simplify the Hamiltonian Hn by replacing the hopping term Jaa (t) by 
its time average. Setting Jaa(t) = Jaa with Jaa = {^/tq) Jq'^ ^tJaait) and rescaling 
the time as t ^ t' — {UinitTQ / AUaa) yields the transformation 

which turns Hf( into the Landau-Zener form Hn = At' + B, where A and B are 
Hcrmitian matrices and A is diagonal [40l [4TJ |42]. The energy spectrum of the 
Hamiltonian (j29p reproduces approximately the features of the spectrum of HcS except 
that the gap is overestimated in the small lattice limit. 

Starting at i = in the small lattice limit where the system is impulse, we use the 
AI approximation to derive the density of defects at the end of a linear ramp which 
drives the system to the large lattice limit, where the system is adiabatic. Defining 
the relaxation time as the inverse of the energy gap 1/AE between the levels of Hf{, 
with AE = [{AUaat/TQ)'^ + (2V2 Jaa)^] ^ , cquation ([23]) can be solved analytically and 
the instant t at which our system exits the impulse regime is given by 



where -q = AJ^^^/ AUaa- We evaluate the density of defects ^ using equation ([M)) with 
l^gi^i)) = Iff(O)) and l^e(i)) = |e(i)). In order to simplify the expression of ^, we 
have set 6(0) = 7r/2, which is a greater value than the one we would obtain using the 
real parameters of the system. This has no significant physical consequences in our 
case as it is actually equivalent to considering a smaller energy gap in the small lattice 
limit. The density of defects is then given by 

f^=i(l-sin^), (31) 
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Figure 9. Comparison between the analytical formula for the density of defects 
(solid line) and the results of the TEBD calculations for the density of kinks 
(dots) 3! = (k)/M, for a system of M = 12 sites with (a) g = g2 and a = ^2 
and (b) g = gi and a = 154. Using the parameters shown in figure |4] we have 

t^mit/S2 = 6.8, AUaa/g2 = 16.4, Jaa/Ea = 0.027. 



where 9 = 9{t). Inverting equation (|31|1 . we obtain 

r,m-'-^4=^, (32) 

which gives an approximate analytical expression of the quench time tq required to 
obtain a density of defects ^. 

In our system, defects correspond mainly to doubly occupied sites. Thus, the 
number of defects is approximately measured by the operator 

M 

k = Y,f^{f^-\). (33) 

1=1 

Hence, the density of defects is given by = {K)/M = ((nf)^) — (nf). Numerical 
calculations show that the density of defects ^ has the same scaling behaviour as 
An^ j . A comparison between equation (|3ip and the numerical values of the density 
of defects in the large lattice limit for different values of the ramping time tq is shown 
in figure [HI For M = 12 particles and g = g2, we find that the analytical formula for 
^ fits the numerical data well for a ~ 1.41. For g ^ gi the fit is less accurate. This 
is expected since for this value of g, the system has a less distinct separation between 
the adiabatic and impulse regime than for g = g2- For the number of particles we 
have been able to simulate, the fit improves as we increase the number of particles 
for both values of g. In addition to this, we see from equation that the time 
required to initialize a register with a given fidelity — and thus the adiabatic time for 
small & — scales with the ratio AUaa/Jia- 



4. Conclusion 



We have shown that the dynamics of an optical lattice whose periodicity is doubled via 
superlattice potentials is very well described by a two-mode Hubbard-like Hamiltonian. 
The parameters of this Hamiltonian have been evaluated in the tight binding 
approximation using optimally localized GWFs. The doubling of the period removes 
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half of the lattice sites and doubles the filling factor. We have shown that this doubling 
can be used for the fast initialization of a quantum register. By starting from a half 
filled lattice in the small lattice limit filled by either a TG {g = gi) gas or a superfiuid 
(ff = 92), a commensurate MI state corresponding to an atomic quantum register can 
obtained on timescales shorter than those achieved by direct quantum freezing of a 
superfiuid with same lattice spacing. Furthermore, we derived an analytical expression 
for the density of defects as a function of the quench time for linear ramping of the 
superlattice. We found that the time required to achieve a given density of defects is 
proportional to the ratio AUaa/Jaa- 

Our numerical calculations of ground state properties suggest that doubling 
the lattice period drives the system through a quantum phase transition for large 
lattices M — > 00. The eventual abrupt change in the ground state properties might 
be observable by time-of-fiight measurements. An investigation of whether such a 
quantum phase transition indeed exists is beyond the scope of the current work but 
will be investigated in future work. 

In this work we concentrated on the transition from filling factor n ~ 1/2 to 
n — 1. We finally note that the idea developed in this paper may be extended by 
considering lattices with an initial filling factor of n = 1/2^ (where i is an integer). 
Subsequently removing every second barrier will create a lattice with period 2a and 
filling factor 1/2^^^. This procedure could be repeated £ times providing a lattice with 
filling factor n = 1 and large lattice spacing 2^a. 
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Appendix 

Appendix A. Definition and localization properties of GWFs 

In the TB limit, the effective single particle Hamiltonian of our system in momentum 
space (in the basis \a)q = [(1/-\/M) J2iLi e''"*dj]|0) with a = a,b) can be written as 



For 11 = and /i = 1, the elements of the matrices £(/i) correspond to the local site 
energy and hopping matrix elements between neighbouring sites, respectively. For the 
TB approximation to be accurate, we need the eigenvalues Ea.q of Hq ^ to reproduce 
very closely the band structure of the exact single particle Hamiltonian Hq^x for all 
values of the lattice profile s. If we were to use WFs as mode functions Wa.i, the 
matrices e would be diagonal [M] and the dispersion relations of the two Bloch bands 




(A.l) 



The elements of the matrices e(/i) are given by [43] 

SapilJ-) = {Wa,i\ho.x\'^l3,i+f,)- 



(A.2) 
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sinusoidal. In order to obtain a more accurate description, we use GWFs as mode 
functions. The definition of GWFs is given by 

(A.3) 

where '^a,q = X]/3=a b ^p'a' ^ P,<! Called the generalized- Bloch orbital with ^a,q the 
Bloch function associated with the band a and Ri is the center of site i [lH [25] . 
The rows of the 2x2 matrix [/^'^^ contain the (real) normalized eigenvectors of Ho,q 
associated with the eigenvalues Ea^q, that is J2fj.=a bi^o,q)n,p,Ua]j, = Ea.qUi% [43] . 
Inserting GWFs in equation (|A.2p . we recover the elements SapifJ.) and, hence, GWFs 
correspond to the mode functions associated with the effective Hamiltonian Ho g |43j . 
Notice that the definition of GWFs reduces to that of WFs for = 1. 

Localization properties of GWFs 

Given a valid set of GWFs, another equally valid set of GWFs can be obtained by 
applying the following transformation on the C/'^'^^ matrices 

- ( "T e*... ) (*^^) 

where 4>a{q) are (real) functions of q which can be chosen freely as long as they 
do not introduce discontinuities in the generalized Bloch function [IS]. The gauge 
transformation (|A.4[) is equivalent to re-phasing each Bloch function as ^'a,g — *■ 
ei<^a('3)\j/^^^. Notice that gauge transformations do not affect the value of the 
parameters Va and Ja,p calculated using the relations ^ and ©, respectively, 
but they alter the localization properties — the spread — of the GWFs. Following the 
convention suggested by Blount |45j . we set the phase functions 4>a{q) in a manner 
that leads to maximally localized GWFs. That is, we choose the phase functions such 
that the resulting GWFs minimize the spread functional 

a 

where in our case a — a,b while (x^) ~ {wa,i\x^\wa,i) and (x)^ = {wa,i\x\wa^i) 
correspond to the center of a GWF and its second moment, respectively. 
Expanding Ua.qix, s) = e~"'^5'Q.q(x, s) into plane waves yields 

Ua,g{x,s)=Y, Ga,j{q,s)e'^^'', (A.6) 
j 

where Kj = 2nj/L with L ~ Ma. Invoking the translational symmetry of the 
lattice and the convolution theorem [45], the value of the functional 17 is minimized 
when the expansion coefficients Ga,j{q-,s) are chosen real, which is always possible 
when the lattice possesses mirror symmetry [43]. This is equivalent to the choice of 
purely real GWFs for even gencralized-Bloch functions and purely imaginary ones for 
odd generalized-Bloch functions. Notice that this conclusion is in agreement with a 
conjecture of Marzari et al. [25] on the real nature (up to a global phase) of maximally 
localized WFs. 

We have numerically evaluated the phases 4)a{Q) using the algorithm described 
in [25] for the special case of ID WFs. This procedure minimized the functional 12 in 
the limit of very fine sampling of the q-space. The effect of this localization procedure 
is illustrated in figure [3^. 
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Appendix B. Parameters of the single-particle Hamiltonian 

We diagonalize Hq in momentum space for periodic boundary conditions. Using the 
Fourier transformations a, = {l/VM)J2g e'«"'aq and U = (l/VM)Eg e'9"*6q, Hq 
becomes block diagonal 

9 

where q ~ 2TTv/Ma, v ~ \ . . . M . In the basis {|a)q, \h)q]. the operator i?o,g reads 

jT Va - 2 Jaa COS qa -i2 Jab sinqa \ , , 

" i2JafcSinga Vb + 2 Jbb cos qa ^ ' 

with \a)q — at|0) and \b)q ~ ^J|0). For simplicity, the explicit dependence of the 
parameters on the lattice profile has been dropped. 

Due to the periodicity of the lattice, the eigenvalues of ifo,g exhibit a band 
structure. By choosing the points q ~ 0, vr/a, 7r/2a in the Brillouin zone, we derive and 
solve a set of equations for the parameters Va and Ja.p as functions of the eigenvalues 

Ea,q of Ho^q 

Ea,TT ~ £-0,0 J _ Eb.O — Eb^-n 



4 



Ea.O ~ Ea.TT Ebfi — Eb,TT Q\ 

Jab=l[{Eb,^-Ea,^f-{Va~Vb) 
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For the eigenvalues of Hq q to reproduce the band structure of the exact single-particle 
Hamiltonian along the x-direction ho^x, we evaluate the parameters Va and Ja,/3 using 
the eigenvalues Sa.q of /lo,^ obtained via exact numerical diagonalization for the same 
points in the Brillouin zone (see figure [J) . The numerical values of the parameters 
obtained via this procedure are shown in figure |4^. 

The accuracy of the TB approximation can be tested by evaluating the standard 
deviation between the exact and approximated band structure. That is, taking 
Nq different points qt on each band, we define the standard deviation between 
the exact and approximate band structure for a given lattice profile by cr^ = 
{l/2Nq) E!I'i(^£a,g. + ^e^?,)' with Aea,q^ = e^.q^ ~ Ea,,q,. Averaging over 
different lattice profiles Si, wc obtain (cr)^ = [(l/A^) ^^^^ 0-^.]^ = 3.4 x 10~^i?fl 
for a lattice depth of Vq = 10 Ej^. This excellent agreement improves further as we 
increase the value of Vq, and hence fully justifies the tight-binding approximation. 



Appendix C. Dynamical and ground state calculations using the TEBD 
algorithm 

The TEBD algorithm is based on directly manipulating a matrix product 
representation of the many-body wave function. Here, we shall briefly describe the 
key aspects of this algorithm and refer the reader to some of the recent literature 
Sa im Sg for more detail. 
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An arbitrary state of a ID quantum lattice system composed of AI sites can be 
written as 

d d 

1^)" (C.l) 

ii=i ij\/=i 

where Cj-^.-.j^^j is a set of d*^ complex amplitudes and \im) is a basis spanning the local 
d-dimensional Hilbert space of site m. Within time-dependent DMRG the amplitudes 
Cjj...jjj are constructed from a product of tensors 

{x} 

{a} = l 

where {a} = {ai, • • • , Q!j\/„i}, {x} = {xii ' ' ' i Xa/-i} and with F and A tensors 
chosen to be constructed from the set of M — 1 Schmidt decompositions for contiguous 
partitions of the system. Specifically, the elements of Aq"' are taken to be Schmidt 
coefficients of the bipartite splitting after site m in j^/^) — X)a=i '^""^ l^a with 
Schmidt rank Xm- The Schmidt states \U^) and spanning the left {1, • • • , m} 

and right {to + 1, • ■ • , M} subsystems of sites respectively are then specified by the 
corresponding sums remaining in equation (jC.2p . 

The usefulness of this representation is based on the observation that for ID 
systems with a Hamiltonian composed of nearest neighbour terms the groundstate 
and low-lying excited states have Schmidt coefficients Aq"^ which rapidly decay with 
a when arranged in descending order. Consequently, rather than allowing the Schmidt 
ranks Xm to grow to their maximum permissible value a much smaller fixed upper- 
limit X can be imposed truncating the representation while still providing a near unit 
overlap with the exact state 1-0). Fixing the Schmidt ranks results in the number of 
parameters scaling as O(dx^-A^) and so curtails the possible exponential growth with 
M seen for general coefficients Cj^...j^^. 

The matrix product representation also permits the efficient update of the state 
after the action of a unitary operator on any two neighboring lattice sites. This 
proceeds by modifying the F tensors associated to the sites and the A tensor linking 
them and requiring a number of operations which scales as 0(d'*x'^)- The resulting 
tensors are then systematically truncated back to a maximum rank of x- 

Dynamical simulations can be performed by decomposing the time evolution 
operator exjp{—iHSt), for small time step St, into a sequence of pairwise unitaries 
via a Suzuki- Trotter expansion. Given the properties outlined such a calculation is 
likely to be accurate for a practical value of x if both the initial state and the states 
generated by the dynamics remain in the low-energy manifold of the system. To 
determine the appropriate x calculations are repeated with increasing values of x 
until the final result converges and are unaffected by further increases. For practical 
purposes the convergence is usually quantified by the robustness of the expectation 
values calculated. The accuracy of a calculation is also gauged by the sum of the 
discarded Schmidt coefficients at each time step - a quantity which should necessarily 
be small - and the deviation of normalization of the final state from unity which 
indicates the accumulated effect of truncation. 

Finally, initial states are typically taken to be the groundstate of the system 
which are found either by applying the DMRG procedure or, as in this work, by 
simulating imaginary time evolution through the repeated application of exp{—Hdt) 
and subsequent rcnormalization of the state. In our simulations, we have used x = 40. 
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Appendix D. The ramp Sgap 

The ramp Sgap can be evaluated as follows. The energy gap between the ground and 
the first excited state associated with a given lattice profile s is given by gap(.s) = ojoe 
where wofc = (eGff,fc—£cff,o) with ereff, a the a-th instantaneous eigenvalue of ifeff(s)- The 
gap is evaluated numerically via the direct diagonalization of HcS for a small number 
of sites. The function Sgap(i) is defined by the equation dsgap(t)/dt — gap(sgap(t))/-Rr 
where = JJ^'' di gap(sgap(t)) is a normalization constant. 
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